Packages
if (!require("tidyverse")) install.packages("tidyverse")
library("tidyverse") #tidyr, ggplot, dplyr, etc
if (!require("ggpubr")) install.packages("ggpubr")
library("ggpubr")
if (!require("lme4")) install.packages("lme4")
library("lme4") # mixed-effect models
if (!require("lmerTest")) install.packages("lmerTest")
library("lmerTest") # p-values
if (!require("car")) install.packages("car")
library("car") # test for VIF
if (!require("broom.mixed")) install.packages("broom.mixed")
library("broom.mixed") # model export
if (!require("emmeans")) install.packages('emmeans')
library('emmeans')
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library("RColorBrewer") # model export
if (!require("rmdformats")) install.packages("rmdformats")
library("rmdformats") # clean html R markdown formatBackground
This code walks through the statistical models and figures for a study on the rate of change of cutaneous evaporative water loss (CEWL) in response to temperature change in Western Fence Lizards (Sceloporus occidentalis). Data was collected and analyzed by Calvin Davis and Savannah Weaver, and the project was advised by Dr. Emily Taylor at California Polytechnic State University. The sections with an asterisk are the stats and figures presented in the publication.
Load Data
Get the RDS files created in the data wrangling Rmd:
lizard_tmt_info <- read_rds("Data/collection_dat_formatted.RDS")
temp_time_series <- read_rds("Data/temp_time_series.RDS") %>%
dplyr::filter(lizard_ID != 410)
CEWL_time_series <- read_rds("Data/CEWL_time_series.RDS") # not needed???
unique(CEWL_time_series$date)## [1] "2021-10-05" "2021-10-11" "2021-10-12" "2021-10-18" "2021-10-19"by_minute_indiv <- temp_time_series %>%
mutate(time_min_block = floor(time_min)) %>%
group_by(date, lizard_ID, treatment, time_min_block) %>%
summarise(CEWL_per_min = mean(CEWL_g_m2h),
clo_temp_per_min = mean(calibrated_cloacal_temp),
dors_temp_per_min = mean(calibrated_dorsal_temp)) ## `summarise()` has grouped output by 'date', 'lizard_ID', 'treatment'. You can
## override using the `.groups` argument.by_minute_avg <- by_minute_indiv %>%
group_by(treatment, time_min_block) %>%
summarise(CEWL_min_avg = mean(CEWL_per_min),
clo_temp_min_avg = mean(clo_temp_per_min),
dors_temp_min_avg = mean(dors_temp_per_min))## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.by_temp_minute_avg <- by_minute_indiv %>%
mutate(clo_temp_chunk = floor(clo_temp_per_min)) %>%
group_by(treatment, clo_temp_chunk) %>%
summarise(CEWL_min_avg = mean(CEWL_per_min))## `summarise()` has grouped output by 'treatment'. You can override using the
## `.groups` argument.Stats & Models
* Temp Correlation
Check the correlation between dorsal and cloacal temps:
##
## Call:
## lm(formula = dors_temp_min_avg ~ clo_temp_min_avg, data = by_minute_avg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3083 -0.5866 0.1572 0.2957 1.6472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.69880 0.65881 10.17 5.19e-13 ***
## clo_temp_min_avg 0.75694 0.02499 30.29 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7538 on 43 degrees of freedom
## Multiple R-squared: 0.9552, Adjusted R-squared: 0.9542
## F-statistic: 917.7 on 1 and 43 DF, p-value: < 2.2e-16CEWL (raw) ~ temp LMM
First, test linear effects of temperature:
CEWL_LMM_02a <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
clo_temp_per_min *
treatment +
(1|date/lizard_ID))
CEWL_LMM_02b <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
dors_temp_per_min *
treatment +
(1|date/lizard_ID))Compare to the inclusion of polynomial factors:
CEWL_LMM_03a <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
poly(clo_temp_per_min, 2)*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_03a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_02a: CEWL_per_min ~ clo_temp_per_min * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_per_min ~ poly(clo_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02a 9 1786.0 1823.1 -883.98 1768.0
## CEWL_LMM_03a 12 1690.4 1740.0 -833.20 1666.4 101.56 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_03b <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
poly(dors_temp_per_min, 2)*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_03b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_02b: CEWL_per_min ~ dors_temp_per_min * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_per_min ~ poly(dors_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_02b 9 1809.6 1846.8 -895.80 1791.6
## CEWL_LMM_03b 12 1724.3 1773.8 -850.15 1700.3 91.304 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The poly models have significantly better explanatory power than the linear models.
Also try log transformation to make sure a polynomial effect is best. But, cannot compare a log-log model because then the response variable would be different.
CEWL_LMM_04a <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
log(clo_temp_per_min)*treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_04a, CEWL_LMM_02a)## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_04a: CEWL_per_min ~ log(clo_temp_per_min) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_02a: CEWL_per_min ~ clo_temp_per_min * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 1801.4 1838.6 -891.72 1783.4
## CEWL_LMM_02a 9 1786.0 1823.1 -883.98 1768.0 15.475 0## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_04a: CEWL_per_min ~ log(clo_temp_per_min) * treatment + (1 | date/lizard_ID)
## CEWL_LMM_03a: CEWL_per_min ~ poly(clo_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04a 9 1801.4 1838.6 -891.72 1783.4
## CEWL_LMM_03a 12 1690.4 1740.0 -833.20 1666.4 117.03 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1CEWL_LMM_04b <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
log(dors_temp_per_min)*treatment +
treatment +
(1|date/lizard_ID))
anova(CEWL_LMM_04b, CEWL_LMM_02b)## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_04b: CEWL_per_min ~ log(dors_temp_per_min) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_02b: CEWL_per_min ~ dors_temp_per_min * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 1825.5 1862.6 -903.73 1807.5
## CEWL_LMM_02b 9 1809.6 1846.8 -895.80 1791.6 15.858 0## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## CEWL_LMM_04b: CEWL_per_min ~ log(dors_temp_per_min) * treatment + treatment + (1 | date/lizard_ID)
## CEWL_LMM_03b: CEWL_per_min ~ poly(dors_temp_per_min, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## CEWL_LMM_04b 9 1825.5 1862.6 -903.73 1807.5
## CEWL_LMM_03b 12 1724.3 1773.8 -850.15 1700.3 107.16 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The logarithmic models are NOT better than the linear or polynomial ones. The poly models are best and the linear and log models are equally less-good.
Check assumptions:
## GVIF Df GVIF^(1/(2*Df))
## poly(clo_temp_per_min, 2) 788782.98078 2 29.801586
## treatment 5.82529 2 1.553565
## poly(clo_temp_per_min, 2):treatment 787225.05195 4 5.457734## GVIF Df GVIF^(1/(2*Df))
## poly(dors_temp_per_min, 2) 4.192428e+05 2 25.445817
## treatment 5.645315e+00 2 1.541424
## poly(dors_temp_per_min, 2):treatment 4.183579e+05 4 5.043053There’s no fanning or different pattern. VIF is negligible.
Re-run with p-values:
CEWL_LMM_besta <- lmerTest::lmer(data = by_minute_indiv,
CEWL_per_min ~
treatment * poly(clo_temp_per_min, 2) +
(1|date/lizard_ID))
summary(CEWL_LMM_besta)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## CEWL_per_min ~ treatment * poly(clo_temp_per_min, 2) + (1 | date/lizard_ID)
## Data: by_minute_indiv
##
## REML criterion at convergence: 1625.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3166 -0.5957 -0.0968 0.5423 5.8100
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 12.210 3.494
## date (Intercept) 24.154 4.915
## Residual 1.578 1.256
## Number of obs: 459, groups: lizard_ID:date, 32; date, 5
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) -0.4704 3.7914 29.3064
## treatmentCooling 8.2628 3.2867 261.5997
## treatmentHeating 17.9752 3.2917 244.1138
## poly(clo_temp_per_min, 2)1 139.9853 50.9856 441.7664
## poly(clo_temp_per_min, 2)2 -502.4452 82.9294 445.3354
## treatmentCooling:poly(clo_temp_per_min, 2)1 -99.1502 51.0336 441.7445
## treatmentHeating:poly(clo_temp_per_min, 2)1 -114.1937 51.0301 441.7380
## treatmentCooling:poly(clo_temp_per_min, 2)2 519.6856 82.9658 445.3323
## treatmentHeating:poly(clo_temp_per_min, 2)2 490.2494 82.9555 445.3375
## t value Pr(>|t|)
## (Intercept) -0.124 0.90211
## treatmentCooling 2.514 0.01254 *
## treatmentHeating 5.461 1.17e-07 ***
## poly(clo_temp_per_min, 2)1 2.746 0.00629 **
## poly(clo_temp_per_min, 2)2 -6.059 2.92e-09 ***
## treatmentCooling:poly(clo_temp_per_min, 2)1 -1.943 0.05267 .
## treatmentHeating:poly(clo_temp_per_min, 2)1 -2.238 0.02573 *
## treatmentCooling:poly(clo_temp_per_min, 2)2 6.264 8.86e-10 ***
## treatmentHeating:poly(clo_temp_per_min, 2)2 5.910 6.82e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnC trtmnH p(___,2)1 p(___,2)2 tC:(___,2)1 tH:(___,2)1
## tretmntClng -0.769
## tretmntHtng -0.764 0.882
## ply(___,2)1 -0.709 0.820 0.814
## ply(___,2)2 0.760 -0.882 -0.872 -0.878
## tC:(___,2)1 0.708 -0.818 -0.814 -0.999 0.877
## tH:(___,2)1 0.708 -0.819 -0.814 -0.999 0.878 0.998
## tC:(___,2)2 -0.760 0.881 0.872 0.878 -1.000 -0.877 -0.877
## tH:(___,2)2 -0.760 0.882 0.872 0.878 -1.000 -0.877 -0.878
## tC:(___,2)2
## tretmntClng
## tretmntHtng
## ply(___,2)1
## ply(___,2)2
## tC:(___,2)1
## tH:(___,2)1
## tC:(___,2)2
## tH:(___,2)2 0.999CEWL_clo_anova <- data.frame(anova(CEWL_LMM_besta, type = "2", ddf = "Kenward-Roger")) %>%
mutate(variable = row.names(.))
CEWL_LMM_bestb <- lmerTest::lmer(data = by_minute_indiv,
CEWL_per_min ~
treatment * poly(dors_temp_per_min, 2) +
(1|date/lizard_ID))
summary(CEWL_LMM_bestb)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: CEWL_per_min ~ treatment * poly(dors_temp_per_min, 2) + (1 |
## date/lizard_ID)
## Data: by_minute_indiv
##
## REML criterion at convergence: 1660.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4453 -0.6089 -0.0814 0.5142 5.9945
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 10.295 3.209
## date (Intercept) 22.052 4.696
## Residual 1.728 1.315
## Number of obs: 459, groups: lizard_ID:date, 32; date, 5
##
## Fixed effects:
## Estimate Std. Error df
## (Intercept) 2.194 3.495 25.854
## treatmentCooling 4.972 2.973 262.431
## treatmentHeating 15.455 2.988 246.451
## poly(dors_temp_per_min, 2)1 179.875 54.756 444.012
## poly(dors_temp_per_min, 2)2 -410.895 72.930 441.411
## treatmentCooling:poly(dors_temp_per_min, 2)1 -149.678 54.791 443.997
## treatmentHeating:poly(dors_temp_per_min, 2)1 -158.197 54.802 443.991
## treatmentCooling:poly(dors_temp_per_min, 2)2 428.187 72.966 441.398
## treatmentHeating:poly(dors_temp_per_min, 2)2 403.389 72.965 441.425
## t value Pr(>|t|)
## (Intercept) 0.628 0.53572
## treatmentCooling 1.672 0.09565 .
## treatmentHeating 5.173 4.78e-07 ***
## poly(dors_temp_per_min, 2)1 3.285 0.00110 **
## poly(dors_temp_per_min, 2)2 -5.634 3.14e-08 ***
## treatmentCooling:poly(dors_temp_per_min, 2)1 -2.732 0.00655 **
## treatmentHeating:poly(dors_temp_per_min, 2)1 -2.887 0.00408 **
## treatmentCooling:poly(dors_temp_per_min, 2)2 5.868 8.65e-09 ***
## treatmentHeating:poly(dors_temp_per_min, 2)2 5.529 5.54e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) trtmnC trtmnH p(___,2)1 p(___,2)2 tC:(___,2)1 tH:(___,2)1
## tretmntClng -0.753
## tretmntHtng -0.748 0.878
## ply(___,2)1 -0.721 0.849 0.843
## ply(___,2)2 0.743 -0.877 -0.868 -0.923
## tC:(___,2)1 0.721 -0.848 -0.843 -0.999 0.923
## tH:(___,2)1 0.720 -0.848 -0.843 -0.999 0.923 0.999
## tC:(___,2)2 -0.743 0.876 0.867 0.923 -1.000 -0.922 -0.922
## tH:(___,2)2 -0.743 0.876 0.867 0.923 -0.999 -0.922 -0.922
## tC:(___,2)2
## tretmntClng
## tretmntHtng
## ply(___,2)1
## ply(___,2)2
## tC:(___,2)1
## tH:(___,2)1
## tC:(___,2)2
## tH:(___,2)2 0.999CEWL_dors_anova <- data.frame(anova(CEWL_LMM_bestb, type = "2", ddf = "Kenward-Roger")) %>%
mutate(variable = row.names(.))
# double check sample sizes
temp_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 11
## 2 Cooling 11
## 3 Heating 10These are the best models for CEWL ~ temperature.
Export:
# model anovas
CEWL_temp_F_stats <- CEWL_clo_anova %>%
rbind(CEWL_dors_anova) %>%
mutate(partial_sum_of_squares = round(Sum.Sq, 0),
F_only = round(F.value, 2),
DenDF = round(DenDF, 0),
F_statistic = paste(F_only, " (", NumDF, ",", DenDF, ")", sep = ""),
p_value = (Pr..F.)) %>%
dplyr::select(variable,
partial_sum_of_squares,
F_statistic,
p_value)
write.csv(CEWL_temp_F_stats,
"./Results_Statistics/CEWL_temp_F_stats.csv")
# model summaries
broom.mixed::tidy(CEWL_LMM_besta) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_clotemp.csv")
broom.mixed::tidy(CEWL_LMM_bestb) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_dorstemp.csv")CEWL (raw) ~ time LMM
NOTE: Start times vary from 97-170 seconds after starting time series.
Make a polynomial model first:
rates_LMM_01 <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
poly(time_min_block, 2) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_01)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
## Data: by_minute_indiv
##
## REML criterion at convergence: 1630.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6540 -0.5335 0.0275 0.4644 5.8089
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.089 3.015
## date (Intercept) 20.204 4.495
## Residual 1.581 1.258
## Number of obs: 459, groups: lizard_ID:date, 32; date, 5
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 15.273 2.211 6.909
## poly(time_min_block, 2)1 -17.822 2.119 -8.409
## poly(time_min_block, 2)2 6.503 2.145 3.032
## treatmentCooling -8.441 1.312 -6.435
## treatmentHeating 2.310 1.361 1.697
## poly(time_min_block, 2)1:treatmentCooling -17.283 3.107 -5.562
## poly(time_min_block, 2)2:treatmentCooling 15.775 3.074 5.132
## poly(time_min_block, 2)1:treatmentHeating 41.968 3.080 13.628
## poly(time_min_block, 2)2:treatmentHeating -14.173 3.082 -4.598
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(t__,2)1 -0.001
## ply(t__,2)2 0.000 0.001
## tretmntClng -0.292 0.002 0.000
## tretmntHtng -0.286 0.002 0.000 0.466
## pl(__,2)1:C 0.001 -0.682 0.000 0.005 -0.001
## pl(__,2)2:C 0.000 0.000 -0.698 0.002 0.000 0.025
## pl(__,2)1:H 0.001 -0.688 0.000 -0.002 -0.001 0.469 0.000
## pl(__,2)2:H 0.000 0.000 -0.696 0.000 0.000 0.000 0.486
## p(__,2)1:H
## ply(t__,2)1
## ply(t__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.008Also make a linear model version:
rates_LMM_02 <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
time_min_block * treatment +
(1|date/lizard_ID))
anova(rates_LMM_02, rates_LMM_01)## refitting model(s) with ML (instead of REML)## Data: by_minute_indiv
## Models:
## rates_LMM_02: CEWL_per_min ~ time_min_block * treatment + (1 | date/lizard_ID)
## rates_LMM_01: CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## rates_LMM_02 9 1786.6 1823.8 -884.31 1768.6
## rates_LMM_01 12 1682.6 1732.2 -829.30 1658.6 110.02 3 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Polynomial is much better than linear.
Also compare to log version:
rates_LMM_03 <- lme4::lmer(data = by_minute_indiv,
CEWL_per_min ~
log(time_min_block) * treatment +
(1|date/lizard_ID))
rates_LMM_03_loglog <- lme4::lmer(data = by_minute_indiv,
log(CEWL_per_min) ~
log(time_min_block) * treatment +
(1|date/lizard_ID))The polynomial model is best (slightly better than log and much better than lm options).
Check linear regression assumptions:
## GVIF Df GVIF^(1/(2*Df))
## poly(time_min_block, 2) 8.011443 2 1.682394
## treatment 1.000098 2 1.000024
## poly(time_min_block, 2):treatment 8.011931 4 1.297081re-run with p-values:
rates_LMM_best <- lmerTest::lmer(data = by_minute_indiv,
CEWL_per_min ~
poly(time_min_block, 2) * treatment +
(1|date/lizard_ID))
summary(rates_LMM_best)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula:
## CEWL_per_min ~ poly(time_min_block, 2) * treatment + (1 | date/lizard_ID)
## Data: by_minute_indiv
##
## REML criterion at convergence: 1630.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6540 -0.5335 0.0275 0.4644 5.8089
##
## Random effects:
## Groups Name Variance Std.Dev.
## lizard_ID:date (Intercept) 9.089 3.015
## date (Intercept) 20.204 4.495
## Residual 1.581 1.258
## Number of obs: 459, groups: lizard_ID:date, 32; date, 5
##
## Fixed effects:
## Estimate Std. Error df t value
## (Intercept) 15.273 2.211 5.075 6.909
## poly(time_min_block, 2)1 -17.822 2.119 421.033 -8.409
## poly(time_min_block, 2)2 6.503 2.145 421.043 3.032
## treatmentCooling -8.441 1.312 25.139 -6.435
## treatmentHeating 2.310 1.361 25.208 1.697
## poly(time_min_block, 2)1:treatmentCooling -17.283 3.107 421.456 -5.562
## poly(time_min_block, 2)2:treatmentCooling 15.775 3.074 421.092 5.132
## poly(time_min_block, 2)1:treatmentHeating 41.968 3.080 421.102 13.628
## poly(time_min_block, 2)2:treatmentHeating -14.173 3.082 421.059 -4.598
## Pr(>|t|)
## (Intercept) 0.000916 ***
## poly(time_min_block, 2)1 6.47e-16 ***
## poly(time_min_block, 2)2 0.002582 **
## treatmentCooling 9.49e-07 ***
## treatmentHeating 0.101948
## poly(time_min_block, 2)1:treatmentCooling 4.74e-08 ***
## poly(time_min_block, 2)2:treatmentCooling 4.38e-07 ***
## poly(time_min_block, 2)1:treatmentHeating < 2e-16 ***
## poly(time_min_block, 2)2:treatmentHeating 5.64e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) pl(__,2)1 pl(__,2)2 trtmnC trtmnH p(__,2)1:C p(__,2)2:C
## ply(t__,2)1 -0.001
## ply(t__,2)2 0.000 0.001
## tretmntClng -0.292 0.002 0.000
## tretmntHtng -0.286 0.002 0.000 0.466
## pl(__,2)1:C 0.001 -0.682 0.000 0.005 -0.001
## pl(__,2)2:C 0.000 0.000 -0.698 0.002 0.000 0.025
## pl(__,2)1:H 0.001 -0.688 0.000 -0.002 -0.001 0.469 0.000
## pl(__,2)2:H 0.000 0.000 -0.696 0.000 0.000 0.000 0.486
## p(__,2)1:H
## ply(t__,2)1
## ply(t__,2)2
## tretmntClng
## tretmntHtng
## pl(__,2)1:C
## pl(__,2)2:C
## pl(__,2)1:H
## pl(__,2)2:H -0.008## Type II Analysis of Variance Table with Kenward-Roger's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## poly(time_min_block, 2) 160.58 80.289 2 421.20 50.768 < 2.2e-16
## treatment 108.67 54.337 2 25.21 34.358 6.314e-08
## poly(time_min_block, 2):treatment 736.17 184.044 4 421.19 116.374 < 2.2e-16
##
## poly(time_min_block, 2) ***
## treatment ***
## poly(time_min_block, 2):treatment ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# double check sample sizes
CEWL_time_series %>%
group_by(lizard_ID, treatment) %>%
summarise(n()) %>%
group_by(treatment) %>%
summarise(n())## `summarise()` has grouped output by 'lizard_ID'. You can override using the
## `.groups` argument.## # A tibble: 3 × 2
## treatment `n()`
## <fct> <int>
## 1 Control 12
## 2 Cooling 12
## 3 Heating 13Export:
broom.mixed::tidy(rates_LMM_best) %>%
mutate(estimate = round(estimate, digits = 2),
std.error = round(std.error, digits = 2),
statistic = round(statistic, digits = 2),
df = round(df, digits = 0)) %>%
dplyr::select(effect, group, term, estimate,
std.error, statistic, df, p.value) %>%
write.csv("./Results_Statistics/CEWL_best_model_time.csv")Rate of Change
* Tb ~ Time
lizard_temp_time_lm_log <- lm(data = temp_time_series,
log(calibrated_cloacal_temp) ~ log(time_min)*lizard_ID)
plot(lizard_temp_time_lm_log)
Now get the trend/slope for each individual and make into a variable in df.
lizard_coeffs_temp <- data.frame(emtrends(lizard_temp_time_lm_log,
"lizard_ID", var = "time_min")) %>%
left_join(lizard_tmt_info,
by = c("lizard_ID" = "individual_ID")) %>%
dplyr::select(lizard_ID, rate_change = time_min.trend,
SE, df, lower.CL, upper.CL,
date,
treatment
) %>%
dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs_temp)## lizard_ID rate_change SE df
## 401 : 1 Min. :-0.0545701 Min. :0.0001623 Min. :24381
## 402 : 1 1st Qu.:-0.0324598 1st Qu.:0.0001752 1st Qu.:24381
## 404 : 1 Median : 0.0021493 Median :0.0001779 Median :24381
## 406 : 1 Mean :-0.0001463 Mean :0.0001863 Mean :24381
## 407 : 1 3rd Qu.: 0.0346028 3rd Qu.:0.0001876 3rd Qu.:24381
## 408 : 1 Max. : 0.0440239 Max. :0.0003102 Max. :24381
## (Other):26
## lower.CL upper.CL date treatment
## Min. :-0.0549235 Min. :-0.0542167 Min. :2021-10-05 Control:11
## 1st Qu.:-0.0328089 1st Qu.:-0.0321107 1st Qu.:2021-10-11 Cooling:11
## Median : 0.0018033 Median : 0.0024954 Median :2021-10-12 Heating:10
## Mean :-0.0005115 Mean : 0.0002189 Mean :2021-10-13
## 3rd Qu.: 0.0342589 3rd Qu.: 0.0349467 3rd Qu.:2021-10-18
## Max. : 0.0436800 Max. : 0.0443679 Max. :2021-10-19
##
## date_factor
## 2021-10-05:6
## 2021-10-11:6
## 2021-10-12:7
## 2021-10-18:7
## 2021-10-19:6
##
## ## [1] 32then get average rate of change by treatment group:
##
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_temp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0170679 -0.0020104 0.0005646 0.0021987 0.0134091
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002336 0.001653 1.413 0.168
## treatmentCooling -0.039838 0.002338 -17.042 < 2e-16 ***
## treatmentHeating 0.035878 0.002395 14.978 3.49e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005482 on 29 degrees of freedom
## Multiple R-squared: 0.9719, Adjusted R-squared: 0.97
## F-statistic: 501.3 on 2 and 29 DF, p-value: < 2.2e-16## treatment emmean SE df lower.CL upper.CL
## 1 Control 0.002336053 0.001652941 29 -0.00104459 0.005716696
## 2 Cooling -0.037502169 0.001652941 29 -0.04088281 -0.034121526
## 3 Heating 0.038214475 0.001733619 29 0.03466883 0.041760123avg_temp_change_pairs <- data.frame(pairs(emmeans(avg_temp_change, "treatment")))
avg_temp_change_pairs## contrast estimate SE df t.ratio p.value
## 1 Control - Cooling 0.03983822 0.002337611 29 17.04228 3.885781e-15
## 2 Control - Heating -0.03587842 0.002395338 29 -14.97844 1.454392e-14
## 3 Cooling - Heating -0.07571664 0.002395338 29 -31.61000 3.552714e-15and finally, t-tests of rate of change per tmt group:
lizard_coeffs_temp_control <- lizard_coeffs_temp %>%
dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_temp_control$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_temp_control$rate_change
## t = 6.5738, df = 10, p-value = 6.282e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.001544266 0.003127840
## sample estimates:
## mean of x
## 0.002336053lizard_coeffs_temp_cool <- lizard_coeffs_temp %>%
dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_temp_cool$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_temp_cool$rate_change
## t = -14.558, df = 10, p-value = 4.66e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.04324210 -0.03176224
## sample estimates:
## mean of x
## -0.03750217lizard_coeffs_temp_heat <- lizard_coeffs_temp %>%
dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_temp_heat$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_temp_heat$rate_change
## t = 32.083, df = 9, p-value = 1.364e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.03552000 0.04090895
## sample estimates:
## mean of x
## 0.03821447* CEWL ~ Time
We want to know why each lizard had a different rate of change of CEWL during the 15-minute temperature treatment and measurement period.
First, calculate rate of change for each lizard:
lizard_CEWL_time_lm_log <- lm(data = temp_time_series,
log(CEWL_g_m2h) ~
log(time_min)*lizard_ID)
plot(lizard_CEWL_time_lm_log)
Now get the trend/slope for each individual and make into a variable in df.
lizard_coeffs_CEWL <- data.frame(emtrends(lizard_CEWL_time_lm_log,
"lizard_ID", var = "time_min")) %>%
left_join(lizard_tmt_info,
by = c("lizard_ID" = "individual_ID")) %>%
dplyr::select(lizard_ID, rate_change = time_min.trend,
SE, df, lower.CL, upper.CL,
date,
treatment,
mass_g, SVL_mm,
chamber_time_elapsed_hr
) %>%
dplyr::mutate(date_factor = as.factor(date))
summary(lizard_coeffs_CEWL)## lizard_ID rate_change SE df
## 401 : 1 Min. :-0.094984 Min. :0.0005118 Min. :24381
## 402 : 1 1st Qu.:-0.033908 1st Qu.:0.0005524 1st Qu.:24381
## 404 : 1 Median :-0.013033 Median :0.0005611 Median :24381
## 406 : 1 Mean :-0.014469 Mean :0.0005875 Mean :24381
## 407 : 1 3rd Qu.: 0.004875 3rd Qu.:0.0005916 3rd Qu.:24381
## 408 : 1 Max. : 0.058490 Max. :0.0009780 Max. :24381
## (Other):26
## lower.CL upper.CL date treatment
## Min. :-0.096089 Min. :-0.093879 Min. :2021-10-05 Control:11
## 1st Qu.:-0.035009 1st Qu.:-0.032807 1st Qu.:2021-10-11 Cooling:11
## Median :-0.014197 Median :-0.011868 Median :2021-10-12 Heating:10
## Mean :-0.015621 Mean :-0.013318 Mean :2021-10-13
## 3rd Qu.: 0.003788 3rd Qu.: 0.005963 3rd Qu.:2021-10-18
## Max. : 0.057419 Max. : 0.059561 Max. :2021-10-19
##
## mass_g SVL_mm chamber_time_elapsed_hr date_factor
## Min. : 8.80 Min. :60.00 Min. :0.9333 2021-10-05:6
## 1st Qu.:10.38 1st Qu.:64.00 1st Qu.:1.0000 2021-10-11:6
## Median :11.85 Median :66.00 Median :1.0250 2021-10-12:7
## Mean :11.68 Mean :66.47 Mean :1.1068 2021-10-18:7
## 3rd Qu.:12.82 3rd Qu.:69.00 3rd Qu.:1.1792 2021-10-19:6
## Max. :15.30 Max. :75.00 Max. :1.5333
## ## [1] 32NOW we can look at how the RATE of CEWL change (CEWL change per MINUTE) is different among lizards.
change_LMM_1 <- lme4::lmer(data = lizard_coeffs_CEWL,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
mass_g + SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_1)## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## rate_change ~ treatment + mass_g + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> -148.61
## treatment 2 -120.58
## mass_g 1 -150.42
## SVL_mm 1 -150.33
## chamber_time_elapsed_hr 1 -150.41## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 0.0279063 0.0139532 39.5378
## mass_g 1 0.0000263 0.0000263 0.0744
## SVL_mm 1 0.0000813 0.0000813 0.2305
## chamber_time_elapsed_hr 1 0.0000538 0.0000538 0.1524Drop mass first based on low SS and resulting AIC:
change_LMM_2 <- lme4::lmer(data = lizard_coeffs_CEWL,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
SVL_mm +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_2)## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## rate_change ~ treatment + SVL_mm + chamber_time_elapsed_hr +
## (1 | date_factor)
## npar AIC
## <none> -150.42
## treatment 2 -121.58
## SVL_mm 1 -152.32
## chamber_time_elapsed_hr 1 -152.19## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 0.0279308 0.0139654 39.8100
## SVL_mm 1 0.0000111 0.0000111 0.0316
## chamber_time_elapsed_hr 1 0.0000680 0.0000680 0.1940drop SVL:
change_LMM_3 <- lme4::lmer(data = lizard_coeffs_CEWL,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment +
chamber_time_elapsed_hr +
# include date as random effect bc sig diff weather
(1|date_factor))
drop1(change_LMM_3)## Single term deletions
##
## Model:
## rate_change ~ treatment + chamber_time_elapsed_hr + (1 | date_factor)
## npar AIC
## <none> -152.32
## treatment 2 -123.20
## chamber_time_elapsed_hr 1 -154.15## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 0.0279241 0.0139620 41.4001
## chamber_time_elapsed_hr 1 0.0000513 0.0000513 0.1521drop chamber_time_elapsed_hr:
change_LMM_4 <- lme4::lmer(data = lizard_coeffs_CEWL,
rate_change ~
treatment +
(1|date_factor))
drop1(change_LMM_4)## boundary (singular) fit: see help('isSingular')## Single term deletions
##
## Model:
## rate_change ~ treatment + (1 | date_factor)
## npar AIC
## <none> -154.15
## treatment 2 -116.99## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## treatment 2 0.027917 0.013959 42.898The simplest model is our best/final model.
Check assumptions:
Assumptions look good.
Save with p-values and remove the random effect because it did not account for any variation:
change_CEWL_best <- lm(data = lizard_coeffs_CEWL,
rate_change ~ # remember estimate is CEWL PER MINUTE
treatment)
summary(change_CEWL_best)##
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_CEWL)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043019 -0.015382 0.001456 0.013347 0.037537
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.009175 0.006084 -1.508 0.14239
## treatmentCooling -0.042790 0.008605 -4.973 2.73e-05 ***
## treatmentHeating 0.030128 0.008817 3.417 0.00189 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02018 on 29 degrees of freedom
## Multiple R-squared: 0.7057, Adjusted R-squared: 0.6854
## F-statistic: 34.77 on 2 and 29 DF, p-value: 1.981e-08Treatment group explained a significant amount of the variation.
and get average change ~ tmt:
## treatment emmean SE df lower.CL upper.CL
## 1 Control -0.009174826 0.006084346 29 -0.021618711 0.003269059
## 2 Cooling -0.051965300 0.006084346 29 -0.064409185 -0.039521415
## 3 Heating 0.020952679 0.006381316 29 0.007901423 0.034003936avg_CEWL_change_pairs <- data.frame(pairs(emmeans(change_CEWL_best, "treatment"))) # p values
avg_CEWL_change_pairs## contrast estimate SE df t.ratio p.value
## 1 Control - Cooling 0.04279047 0.008604565 29 4.972997 7.921430e-05
## 2 Control - Heating -0.03012750 0.008817055 29 -3.416958 5.223942e-03
## 3 Cooling - Heating -0.07291798 0.008817055 29 -8.270106 1.205048e-08and finally, t-tests of rate of change per tmt group:
lizard_coeffs_CEWL_control <- lizard_coeffs_CEWL %>%
dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_CEWL_control$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_control$rate_change
## t = -2.64, df = 10, p-value = 0.02473
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.016918337 -0.001431314
## sample estimates:
## mean of x
## -0.009174826lizard_coeffs_CEWL_cool <- lizard_coeffs_CEWL %>%
dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_CEWL_cool$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_cool$rate_change
## t = -6.9758, df = 10, p-value = 3.826e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.06856350 -0.03536711
## sample estimates:
## mean of x
## -0.0519653lizard_coeffs_CEWL_heat <- lizard_coeffs_CEWL %>%
dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_CEWL_heat$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_heat$rate_change
## t = 3.0047, df = 9, p-value = 0.01484
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.005178089 0.036727269
## sample estimates:
## mean of x
## 0.02095268* CEWL ~ Temp
lizard_CEWL_temp_loglog_lm <- lm(data = temp_time_series,
log(CEWL_g_m2h) ~ log(calibrated_cloacal_temp) * lizard_ID)
plot(lizard_CEWL_temp_loglog_lm)
Now get the trend/slope for each individual and make into a variable in df.
lizard_coeffs_CEWL_temp <- data.frame(emtrends(lizard_CEWL_temp_loglog_lm,
"lizard_ID", var = "calibrated_cloacal_temp")) %>%
left_join(lizard_tmt_info,
by = c("lizard_ID" = "individual_ID")) %>%
dplyr::select(lizard_ID, rate_change = calibrated_cloacal_temp.trend,
SE, df, lower.CL, upper.CL,
date,
treatment
) %>%
dplyr::mutate(date_factor = as.factor(date))
length(unique(lizard_coeffs_CEWL_temp$lizard_ID))## [1] 32then get average rate of change by treatment group:
avg_CEWL_temp_change <- lm(data = lizard_coeffs_CEWL_temp, rate_change ~ treatment)
summary(avg_CEWL_temp_change)##
## Call:
## lm(formula = rate_change ~ treatment, data = lizard_coeffs_CEWL_temp)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42655 -0.01854 -0.00098 0.03850 0.33076
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.15748 0.04090 -3.850 0.000600 ***
## treatmentCooling 0.21314 0.05784 3.685 0.000935 ***
## treatmentHeating 0.17726 0.05927 2.991 0.005630 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1357 on 29 degrees of freedom
## Multiple R-squared: 0.3475, Adjusted R-squared: 0.3025
## F-statistic: 7.722 on 2 and 29 DF, p-value: 0.002048avg_CEWL_temp_change_coeff <- data.frame(emmeans(avg_CEWL_temp_change, "treatment"))
avg_CEWL_temp_change_coeff## treatment emmean SE df lower.CL upper.CL
## 1 Control -0.15747877 0.04090155 29 -0.24113183 -0.07382571
## 2 Cooling 0.05566413 0.04090155 29 -0.02798893 0.13931719
## 3 Heating 0.01977732 0.04289790 29 -0.06795875 0.10751338avg_CEWL_temp_change_pairs <- data.frame(pairs(emmeans(avg_CEWL_temp_change, "treatment")))
avg_CEWL_temp_change_pairs## contrast estimate SE df t.ratio p.value
## 1 Control - Cooling -0.21314290 0.05784352 29 -3.6848188 0.002611154
## 2 Control - Heating -0.17725609 0.05927197 29 -2.9905549 0.015087163
## 3 Cooling - Heating 0.03588681 0.05927197 29 0.6054601 0.818280857and finally, t-tests of rate of change per tmt group:
lizard_coeffs_CEWL_temp_control <- lizard_coeffs_CEWL_temp %>%
dplyr::filter(treatment == "Control")
t.test(lizard_coeffs_CEWL_temp_control$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_temp_control$rate_change
## t = -2.3009, df = 10, p-value = 0.04419
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -0.309978692 -0.004978853
## sample estimates:
## mean of x
## -0.1574788lizard_coeffs_CEWL_temp_cool <- lizard_coeffs_CEWL_temp %>%
dplyr::filter(treatment == "Cooling")
t.test(lizard_coeffs_CEWL_temp_cool$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_temp_cool$rate_change
## t = 4.905, df = 10, p-value = 0.0006185
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.03037812 0.08095015
## sample estimates:
## mean of x
## 0.05566413lizard_coeffs_CEWL_temp_heat <- lizard_coeffs_CEWL_temp %>%
dplyr::filter(treatment == "Heating")
t.test(lizard_coeffs_CEWL_temp_heat$rate_change)##
## One Sample t-test
##
## data: lizard_coeffs_CEWL_temp_heat$rate_change
## t = 2.8899, df = 9, p-value = 0.01789
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.004295961 0.035258674
## sample estimates:
## mean of x
## 0.01977732Figures
COLOR
cloacal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = calibrated_cloacal_temp,
group = lizard_ID,
color = treatment) +
geom_line(size = 0.2) +
geom_smooth(se=F, method = "lm") +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = "Internal Body Temperature (°C)",
color = "Treatment", se = FALSE) -> ctemp_time## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.## `geom_smooth()` using formula = 'y ~ x'
ggplot(temp_time_series) +
aes(x = log(time_sec),
y = log(calibrated_cloacal_temp),
group = lizard_ID,
color = treatment) +
geom_line() +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = "Internal Body Temperature (°C)",
color = "Treatment", se = FALSE) -> ctemp_time_log
ctemp_time_log
* cloacal temp (relative) ~ time
ggplot(temp_time_series) +
aes(x = time_min,
y = relative_clo_temp,
group = lizard_ID,
color = treatment) +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_line() +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 1, r = 1, b = 1, l = 4.3, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_x_continuous(limits = c(1,15),
breaks = c(seq(1, 15, 2)),
labels = c(seq(1, 15, 2))) +
scale_y_continuous(limits = c(-21,21),
breaks = c(-20, -10, 0, 10, 20),
labels = c(-20, -10, 0, 10, 20)) +
annotate("text", x = 12, y = -2, label = "Control", parse = T,
size = 2.9, family = "sans", color = control) +
annotate("text", x = 4, y = -11, label = "Cooling", parse = T,
size = 2.9, family = "sans", color = cool) +
annotate("text", x = 4, y = 10, label = "Heating", parse = T,
size = 2.9, family = "sans", color = heat) +
labs(x = "Time (minutes)" ,
y = expression('Relative T'[b]~' (°C)'),
color = "Treatment", se = FALSE) -> ctemp_time_relative
ctemp_time_relative
* Temp Rate of Change
ggplot() +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_point(data = lizard_coeffs_temp,
aes(x = treatment,
y = rate_change,
shape = treatment,
color = treatment,
fill = treatment),
size = 1,
alpha = 0.3,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = avg_temp_change_coeff,
aes(x = treatment,
y = emmean,
color = treatment,
group = treatment,
ymin = lower.CL,
ymax = upper.CL),
width = 0.3,
alpha = 0.9) +
geom_point(data = avg_temp_change_coeff,
aes(x = treatment,
y = emmean,
shape = treatment,
fill = treatment),
color = "black",
size = 3) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 1, r = 1, b = 1, l = 8, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
ylim(-1.8, 1.8) +
annotate("text", x = 1, y = 1.8, label = "a", size = 3) +
annotate("text", x = 2, y = 1.8, label = "b", size = 3) +
annotate("text", x = 3, y = 1.8, label = "c", size = 3) +
annotate("text", x = 1, y = 1.55, label = "***", size = 4) +
annotate("text", x = 2, y = 1.55, label = "***", size = 4) +
annotate("text", x = 3, y = 1.55, label = "***", size = 4) +
labs(x = "" ,
y = expression(Delta ~ 'log-T'[b]~' log-'*min^-1*''),
color = "Treatment") -> rate_change_temp_fig
rate_change_temp_fig
dorsal temp ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = calibrated_dorsal_temp,
group = lizard_ID,
color = treatment) +
geom_line() +
theme_classic() +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x= "Time (seconds)" ,
y= "Surface Body Temperature (°C)",
color ="Treatment") -> dtemp_time
dtemp_time
CEWL (raw) ~ time
ggplot(temp_time_series) +
aes(x = time_sec,
y = CEWL_g_m2h,
group = lizard_ID,
color = treatment) +
geom_line() +
# geom_point(size = 0.01,
# alpha = 0.25) +
# geom_smooth(method = "loess",
# size = 2.5,
# se = F) +
theme_classic() +
theme(text = element_text(size = 15)) +
theme(axis.text.x = element_text(size = 10)) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = bquote('CEWL (g/'*m^2*'h)'),
color = "Treatment",
se = FALSE) +
scale_x_continuous(limits = c(0, 925),
breaks=seq(0, 950, 100)) -> CEWL_time
CEWL_time
CEWL (MIN) ~ time
ggplot(by_minute_avg) +
aes(x= time_min_block,
y = CEWL_min_avg,
group = treatment,
color = treatment
) +
geom_point(size = 1,
alpha = 0.5) +
geom_line() +
geom_smooth(method = "lm",
se = F) +
theme_classic() +
theme(text = element_text(size = 15)) +
theme(axis.text.x = element_text(size = 10)) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time (seconds)" ,
y = bquote('CEWL (g/'*m^2*'h)'),
color = "Treatment",
se = FALSE) -> CEWL_time_min
CEWL_time_min## `geom_smooth()` using formula = 'y ~ x'
legend
This code uses one of our plots to make a large, detailed legend, which we can save and use in ggarrange later.
ggplot(temp_time_series) +
geom_smooth(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
color = treatment),
method = "lm",
formula = y ~ poly(x, 2),
size = 1,
se = F) +
geom_point(aes(x = (calibrated_cloacal_temp),
y = CEWL_g_m2h,
fill = treatment,
shape = treatment),
size = 4,
alpha = 1) +
scale_x_continuous(limits = c(14, 39),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = c(0.5,0.5)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") -> fancy_legend
fancy_legend # view
* CEWL (raw) ~ cloacal temperature
ggplot(temp_time_series) +
aes(
x = calibrated_cloacal_temp,
y = CEWL_g_m2h,
group = lizard_ID,
color = treatment
) +
geom_line() +
# annotations for cooling group
annotate(geom = "point", x = 13, y = 3,
size = 3, pch = 25, fill = cool) +
annotate("text", x = 11, y = 2.5, label = "T[15]", parse = T,
size = 3, family = "sans", color = cool) +
annotate(geom = "point", x = 36, y = 19,
size = 3, pch = 25, fill = cool) +
annotate("text", x = 38, y = 18.5, label = "T[1]", parse = T,
size = 3, family = "sans", color = cool) +
# annotations for heating group
annotate(geom = "point", x = 15, y = 14,
size = 3, pch = 24, fill = heat) +
annotate("text", x = 13, y = 14, label = "T[1]", parse = T,
size = 3, family = "sans", color = heat) +
annotate(geom = "point", x = 38, y = 30,
size = 3, pch = 24, fill = heat) +
annotate("text", x = 40, y = 30, label = "T[15]", parse = T,
size = 3, family = "sans", color = heat) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_x_continuous(limits = c(10, 40),
breaks = c(seq(10, 40, 10)),
labels = c(seq(10, 40, 10))) +
scale_y_continuous(limits = c(0, 31),
breaks = c(seq(0, 30, 10)),
labels = c(seq(0, 30, 10))) +
labs(x = expression('T'[b]~' (°C)'),
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE) -> CEWL_ctemp
CEWL_ctemp
CEWL (log-log) ~ cloacal temperature
ggplot(temp_time_series) +
aes(x = log(calibrated_cloacal_temp),
y = log(CEWL_g_m2h),
group = lizard_ID,
shape = treatment,
color = treatment) +
geom_line() +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "right",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
labs(x = "log- Internal Body Temperature (°C)" ,
y = bquote('log- CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE) -> CEWL_ctemp_log
CEWL_ctemp_log
CEWL (MIN) ~ cloacal temperature
ggplot(by_minute_indiv) +
aes(
x = log(clo_temp_per_min),
y = log(CEWL_per_min),
group = lizard_ID,
shape = treatment,
color = treatment
) +
geom_point(size = 1,
alpha = 0.8) +
geom_line() +
# scale_x_continuous(limits = c(12.5, 40),
# breaks = c(seq(15, 35, 5)),
# labels = c(seq(15, 35, 5))) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE)
CEWL (MIN avg) ~ cloacal temperature
ggplot(by_minute_avg) +
aes(
x = clo_temp_min_avg,
y = CEWL_min_avg,
#group = lizard_ID,
shape = treatment,
color = treatment
) +
geom_point(size = 1,
alpha = 0.8) +
geom_line() +
geom_smooth(method = "lm",
#formula = y ~ poly(x, 2),
size = 1.5,
se = F) +
scale_x_continuous(limits = c(12.5, 40),
breaks = c(seq(15, 35, 5)),
labels = c(seq(15, 35, 5))) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
se = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggplot(by_temp_minute_avg) +
aes(
x = clo_temp_chunk,
y = CEWL_min_avg,
#group = lizard_ID,
shape = treatment,
color = treatment
) +
geom_point(size = 1,
alpha = 0.8) +
geom_line() +
geom_smooth(method = "lm",
#formula = y ~ poly(x, 2),
size = 1.5,
se = F)## `geom_smooth()` using formula = 'y ~ x'
CEWL Control
# get control data only
temp_time_series %>%
dplyr::filter(treatment == "Control") %>%
# pipe into ggplot
ggplot() +
aes(x = calibrated_cloacal_temp,
y = CEWL_g_m2h,
color = lizard_ID) +
geom_point(size = 0.1,
alpha = 0.1) +
geom_smooth(method = "lm",
size = 1,
se = F) +
ylim(0,32) +
xlim(25,31) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text.align = 0,
legend.position = "none",
plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt")) +
scale_color_manual(values = cntrls) +
labs(x = "Internal Body Temperature (°C)" ,
y = bquote('CEWL (g '*m^-2*' '*h^-1*')'),
color = "Treatment",
se = FALSE) -> CEWL_control
CEWL_control## `geom_smooth()` using formula = 'y ~ x'
# save
ggsave(filename = "CEWL_temp_control_lm.pdf",
plot = CEWL_control,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 80, height = 60)## `geom_smooth()` using formula = 'y ~ x'* CEWL (relative) ~ time
# check n lizard sample size
temp_time_series %>%
dplyr::filter(complete.cases(time_min, relative_CEWL)) %>%
group_by(lizard_ID) %>%
summarise(max(relative_CEWL)) %>%
nrow()## [1] 32# plot
ggplot(temp_time_series) +
aes(
x= time_min,
y = relative_CEWL,
group = lizard_ID,
color = treatment
) +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_line() +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.title.y = element_text(vjust = 0.5),
plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat),
name = "Treatment") +
scale_x_continuous(limits = c(1,15),
breaks = c(seq(1, 15, 2)),
labels = c(seq(1, 15, 2))) +
scale_y_continuous(limits = c(-18,18),
breaks = c(-10, 0, 10),
labels = c(-10, 0, 10)) +
# annotate("text", x = 12, y = -2, label = "Control", parse = T, # too cluttered
# size = 5, family = "sans", color = control) +
annotate("text", x = 4, y = -15, label = "Cooling", parse = T, #9.3
size = 2.9, family = "sans", color = cool) +
annotate("text", x = 4, y = 7, label = "Heating", parse = T,
size = 2.9, family = "sans", color = heat) +
labs(x = " ", #"Time (minutes)" ,
y = bquote('Relative CEWL (g '*m^-2*' '*h^-1*')')) -> CEWLr_time
CEWLr_time
* CEWL Rate of Change
ggplot() +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_point(data = lizard_coeffs_CEWL,
aes(x = treatment,
y = rate_change,
shape = treatment,
color = treatment,
fill = treatment),
size = 1,
alpha = 0.3,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = avg_CEWL_change_coeff,
aes(x = treatment,
y = emmean,
color = treatment,
group = treatment,
ymin = lower.CL,
ymax = upper.CL),
width = .3,
alpha = 0.9) +
geom_point(data = avg_CEWL_change_coeff,
aes(x = treatment,
y = emmean,
shape = treatment,
fill = treatment),
color = "black",
size = 3) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.title.y = element_text(vjust = -1.8),
plot.margin = margin(t = 1, r = 1, b = 1, l = 1, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_y_continuous(limits = c(-0.1, 0.1),
breaks = c(-0.1, 0, 0.1),
labels = c(-0.1, 0, 0.1)) +
annotate("text", x = 1, y = 0.1, label = "a", size = 3) +
annotate("text", x = 2, y = 0.1, label = "b", size = 3) +
annotate("text", x = 3, y = 0.1, label = "c", size = 3) +
annotate("text", x = 1, y = 0.086, label = "*", size = 4) +
annotate("text", x = 2, y = 0.086, label = "***", size = 4) +
annotate("text", x = 3, y = 0.086, label = "*", size = 4) +
labs(x = "" ,
y = expression(Delta ~ 'log-CEWL log-'*min^-1*''),
color = "Treatment") -> rate_change_CEWL_fig
rate_change_CEWL_fig
* CEWL-Temp Rate of Change
ggplot() +
geom_hline(yintercept = 0, size = 0.5,
linetype = "dashed", color = "black") +
geom_point(data = lizard_coeffs_CEWL_temp,
aes(x = treatment,
y = rate_change,
shape = treatment,
color = treatment,
fill = treatment),
size = 1,
alpha = 0.3,
position = position_jitter(height = 0, width = 0.2)) +
geom_errorbar(data = avg_CEWL_temp_change_coeff,
aes(x = treatment,
y = emmean,
color = treatment,
group = treatment,
ymin = lower.CL,
ymax = upper.CL),
width = .3,
alpha = 0.9) +
geom_point(data = avg_CEWL_temp_change_coeff,
aes(x = treatment,
y = emmean,
shape = treatment,
fill = treatment),
color = "black",
size = 3) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.title.y = element_text(vjust = -0.5),
plot.margin = margin(t = 1, r = 1, b = 4, l = 1, unit = "pt")) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_fill_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
scale_shape_manual(values = c(21,25,24),
name = "Treatment") +
scale_y_continuous(limits = c(-0.6, 0.6),
breaks = c(-0.5, 0, 0.5),
labels = c(-0.5, 0, 0.5)) +
# pairwise comparisons between tmt groups
annotate("text", x = 1, y = 0.6, label = "a", size = 3) +
annotate("text", x = 2, y = 0.6, label = "b", size = 3) +
annotate("text", x = 3, y = 0.6, label = "b", size = 3) +
# t-tests for each tmt group
annotate("text", x = 1, y = 0.52, label = "*", size = 4) +
annotate("text", x = 2, y = 0.52, label = "**", size = 4) +
annotate("text", x = 3, y = 0.52, label = "*", size = 4) +
labs(x = "" ,
y = expression(Delta ~ 'log-CEWL log-T'[b]),
color = "Treatment") -> rate_change_CEWL_temp_fig
rate_change_CEWL_temp_fig
* Experimental Design Skematic (fig 1)
We created reference values to use to make a figure that would show each step of our experiment.
ggplot(exp_design) +
aes(x = time_min,
y = temp_C,
color = treatment) +
geom_line(size = 1) +
geom_vline(xintercept = 60,
#alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(30, 38, family = "sans"), label = "incubation (a)",
color = "black", size = 2.5) +
geom_vline(xintercept = 65,
#alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(62.5, 38, family = "sans"), label = "(b)",
color = "black", size = 2.5) +
geom_vline(xintercept = 80,
#alpha = 0.5,
size = 0.5,
linetype = "dashed", color = "black") +
geom_text(aes(72.5, 38, family = "sans"), label = "tmt (c)",
color = "black", size = 2.5) +
theme_classic() +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 8),
legend.text = element_text(color = "black",
family = "sans",
size = 12),
legend.position = "none",
axis.text.x = element_blank(),
#axis.line.x = element_blank(),
axis.ticks.x = element_blank()) +
scale_color_manual(values = c("Control" = control,
"Cooling" = cool,
"Heating" = heat)) +
labs(x = "Time" ,
y = expression('Target T'[b]~' (°C)'),
color = "Treatment") +
annotate(geom = "point", x = 80, y = 25,
size = 3, pch = 21, fill = control) +
annotate(geom = "point", x = 80, y = 20,
size = 3, pch = 25, fill = cool) +
annotate(geom = "point", x = 80, y = 30,
size = 3, pch = 24, fill = heat) +
annotate("text", x = 35, y = 24, label = "Control", parse = T,
size = 2.9, family = "sans", color = control) +
annotate("text", x = 18, y = 30, label = "Cooling", parse = T,
size = 2.9, family = "sans", color = cool) +
annotate("text", x = 18, y = 20, label = "Heating", parse = T,
size = 2.9, family = "sans", color = heat) +
scale_y_continuous(limits = c(15, 38),
breaks=seq(15, 35, 5)) -> methods_schmatic_fig
methods_schmatic_fig
Save Grouped Figures
* CEWL & Tb x relative & rate of change (Fig 2)
CEWL_temp_multi <- ggarrange(
CEWLr_time, rate_change_CEWL_fig, ctemp_time_relative, rate_change_temp_fig,
nrow = 2, ncol = 2,
widths = c(2,1.2),
labels = c("A", "B", "C", "D"),
label.x = 0.01,
label.y = 1.015
)
ggsave(filename = "CEWL_temp_relative_and_change.pdf",
plot = CEWL_temp_multi,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 160, height = 140)* CEWL ~ Tb
CEWL_temp_relation <- ggarrange(
CEWL_ctemp, rate_change_CEWL_temp_fig,
nrow = 1, ncol = 2,
widths = c(2,1.2),
labels = c("A", "B"),
#label.x = 0.01,
label.y = 1.015
)
ggsave(filename = "CEWL_versus_temp.pdf",
plot = CEWL_temp_relation,
path = "./Results_Figures",
device = "pdf",
dpi = 600,
units = "mm",
width = 160, height = 70)Relative CEWL & Rate of Change
multi_fig_relative <- ggarrange(CEWLr_time,
ggarrange(rate_change_CEWL_fig, LEGEND,
nrow = 1, ncol = 2,
widths = c(2,1),
align = "hv"
),
nrow = 2, ncol = 1,
heights = c(2,1),
labels = c("A", "B"))
multi_fig_relative
